//
// This file is part of cpp-lib, a library of C++ functions and classes.
// Cpp-lib is copyright (C) 2004 and onwards, The CPL Interest Group.
//
// This library is free software; you can redistribute it and/or modify it
// under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation; either version 2.1 of the License, or (at
// your option) any later version.
//
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
// General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with this library (see the file COPYING); if not, write to the
// Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
// 02111-1307  USA
//
// $Id$
//


#ifndef CPP_LIB_MATRIX_H
#define CPP_LIB_MATRIX_H

#include <vector>
#include <functional>
#include <cassert>
#include <iosfwd>
#include <algorithm>
#include <limits>
#include <string>
#include <fstream>

#include <cmath>


/// Expression template based Linear Algebra package, mathematical
/// utilities, etc.

///
/// WARNING.  This package does NO aliasing checks whatsoever.  If you
/// write, for example,
///
///   x = M * x ;
///
/// be aware that while the assignment is in progress, the values on the
/// right hand side will already have changed.  Thus, consider writing
/// instead:
///
///   y = M * x ; x = y ;
///
/// to achieve the desired effect.
///
/// WARNING 2.  Always keep in mind that all functions with matrix
/// return value actually return expression templates with delayed
/// evaluation.  This fact will manifest itself dramatically if you
/// write products of the kind
///
///   M = M_1 * M_2 * ... * M_n
///
/// where n >= 3.
///
/// As above, you should explicitly introduce variables to hold the
/// intermediate result(s).
///

namespace cpl {

namespace math {

/// \defgroup tags Matrix expression tags.
///
//@{

/// \defgroup basic Basic matrix expressions actually containing values.
///
//@{

/// Dense %matrix tag.
struct dense {} ;

/// Dense column vector tag.
struct dense_column {} ;

/// Banded %matrix tag.
struct banded {} ;

/// Fixed %matrix tag.  Dimensions are fixed at compile time.
template< long rows , long columns > struct fixed {} ;

//@}

/// \defgroup proxies Tags for matrix expression templates.
///
//@{

/// Transposed matrix.
template< typename tag > struct transposed {} ;

/// Submatrix defined by row range (assignable).
template< typename tag > struct row_range {} ;

/// Submatrix defined by row range (immutable).
template< typename tag > struct const_row_range {} ;

/// Binary operator op applied to rows * columns index pairs.
template< typename op , long rows , long columns > struct index_fun {} ;

/// Unary operator applied to each element.
template< typename tag , typename op > struct un_op {} ;

/// Matrix product.
template< typename tag1 , typename tag2 > struct product {} ;
/// Row scaling.
template< typename tag1 , typename tag2 > struct row_scaled {} ;
/// Elementwise binary operation with two matrix arguments.
template< typename tag1 , typename tag2 , typename op > struct bin_op {} ;

/// View of a std::vector< T > as a matrix.
struct vector {} ;

/// View of a std::vector< std::vector< T > > as a matrix.
struct vector_vector {} ;

//@} // end proxies

//@}


/// Helper structure for initialize.  Never to be used directly!

template< typename T >
struct initializer {

  initializer( T const& data , long const size )
  : data( data ) , size( size ) {}

  T const& data ;
  long const size ;

} ;

/// \defgroup matrix_templates Templates for various matrix expressions.
///
//@{


/// The general matrix type.

/// T is the type of the contained values, tag is the
/// type of the matrix expression.

template< typename T , typename tag > struct matrix ;


/// Return true if A has entries out of the nonzero band.

template< typename T , typename tag >
inline bool has_out_of_band_entries( matrix< T , tag > const& A ) {

  return    lower_bandwidth( A ) < A.rows() - 1
         || upper_bandwidth( A ) < A.columns() - 1 ;

}


/// Assert() that row index is within bounds.

template< typename T , typename tag >
inline void check_row_index( matrix< T , tag > const& A , long i )
{ assert( 1 <= i && i <= A.rows() ) ; }

/// Assert() that column index is within bounds.

template< typename T , typename tag >
inline void check_column_index( matrix< T , tag > const& A , long j )
{ assert( 1 <= j && j <= A.columns() ) ; }

/// Assert() that row and column index are within bounds.

template< typename T , typename tag >
inline void check_indices( matrix< T , tag > const& A , long i , long j ) {

  check_row_index( A , i ) ;
  check_column_index( A , j ) ;

}

/// Return zero-based index into data array for row major layout.
/// Asserts that \a i, \a j >= 1 and j <= columns.

inline unsigned long row_major_index
( long const i , long const j , long const columns ) ;


//
// Traits structures.
//

// Traits structure defining bandwidth properties of the matrix in
// question.

template< typename T , typename tag > struct bandwidth_traits {

  // Not to be called directly but by the respective free function
  // templates!
  long const lower_bandwidth_( matrix< T , tag > const& A ) const 
  { return A.rows()    - 1 ; }
  long const upper_bandwidth_( matrix< T , tag > const& A ) const 
  { return A.columns() - 1 ; }

} ;


// The canonical way to determine bandwidth is through these function
// templates!

template< typename T , typename tag >
inline long const lower_bandwidth( matrix< T , tag > const& A ) 
{ return bandwidth_traits< T , tag >().lower_bandwidth_( A ) ; }

template< typename T , typename tag >
inline long const upper_bandwidth( matrix< T , tag > const& A ) 
{ return bandwidth_traits< T , tag >().upper_bandwidth_( A ) ; }


template< typename T >
struct bandwidth_traits< T , banded > {

  long const lower_bandwidth_
  ( matrix< T , banded > const& A ) const { return A.get_lower_bandwidth() ; }
  
  long const upper_bandwidth_
  ( matrix< T , banded > const& A ) const { return A.get_upper_bandwidth() ; }

} ;


template< typename T , typename op , long rows , long columns >
struct bandwidth_traits< T , index_fun< op , rows , columns > > {

  long const lower_bandwidth_
  ( matrix< T , index_fun< op , rows , columns > > const& A )
  { return A.function().lower_bandwidth() ; }

  long const upper_bandwidth_
  ( matrix< T , index_fun< op , rows , columns > > const& A )
  { return A.function().upper_bandwidth() ; }

} ;


template< typename T , typename tag > 
struct bandwidth_traits< T , transposed< tag > > {

  long const lower_bandwidth_( matrix< T , transposed< tag > > const& A ) const
  { return upper_bandwidth( A.contained() ) ; }

  long const upper_bandwidth_( matrix< T , transposed< tag > > const& A ) const
  { return lower_bandwidth( A.contained() ) ; }

} ;


template< typename T , typename tag1 , typename tag2 >
struct bandwidth_traits< T , product< tag1 , tag2 > > {

  long const lower_bandwidth_
  ( matrix< T , product< tag1 , tag2 > > const& A ) const {
    return std::min( 
        lower_bandwidth( A.left_operand () ) 
      + lower_bandwidth( A.right_operand() ) , 
      A.rows   () - 1 
    ) ; 
  }
  
  long const upper_bandwidth_
  ( matrix< T , product< tag1 , tag2 > > const& A ) const {
    return std::min( 
        upper_bandwidth( A.left_operand () ) 
      + upper_bandwidth( A.right_operand() ) , 
      A.columns() - 1 
    ) ; 
  }

} ;
 
 

template< typename T , typename tag1 , typename tag2 , typename op >
struct bandwidth_traits< T , bin_op< tag1 , tag2 , op > > {

  long const lower_bandwidth_
  ( matrix< T , bin_op< tag1 , tag2 , op > > const& A ) const
  { return std::max( 
      lower_bandwidth( A.left_operand () ) , 
      lower_bandwidth( A.right_operand() )
    ) ;
  }
  
  long const upper_bandwidth_
  ( matrix< T , bin_op< tag1 , tag2 , op > > const& A ) const
  { return std::max( 
      upper_bandwidth( A.left_operand () ) , 
      upper_bandwidth( A.right_operand() )
    ) ;
  }

} ;


template< typename T , typename tag , typename op >
struct bandwidth_traits< T , un_op< tag , op > > {

  long const lower_bandwidth_
  ( matrix< T , un_op< tag , op > > const& A ) const
  { assert( A.function()( 0 ) == 0 ) ; 
    return lower_bandwidth( A.operand() ) ; }
  
  long const upper_bandwidth_
  ( matrix< T , un_op< tag , op > > const& A ) const
  { assert( A.function()( 0 ) == 0 ) ; 
    return upper_bandwidth( A.operand() ) ; }

} ;


/// Specialization for dense column vectors.

template< typename T >
struct matrix< T , dense_column > {

  typedef T      &       reference ;
  typedef T const& const_reference ;

  matrix( long const columns )
  : data( ( assert( columns >= 1 ) , columns ) )
  {}

  /// Default constructor makes an empty matrix (zero rows and columns).

  matrix() : data() {}

  /// Templated copy constructor delegating to assignment operator.
  /// \todo This is somewhat inefficient.  The members are initialized and
  /// immediately overwritten.

  template< typename tag >
  matrix( matrix< T , tag > const& A )
  : data()
  { *this = A ; }

  /// Swap with other matrix

  void swap( matrix< T , dense >& other ) { data.swap( other.data ) ; }

  /// Set all elements to specified value.

  void fill( T const& x ) { std::fill( data.begin() , data.end() , x ) ; }

  long const rows   () const { return 1           ; }
  long const columns() const { return data.size() ; }

        reference operator()( long const i , long const j = 1 )
  { check_indices( *this , i , j ) ;
    return data[ i - 1 ] ; }

  const_reference operator()( long const i , long const j = 1 ) const
  { check_indices( *this , i , j ) ;
    return data[ i - 1 ] ; }

  /// Templated assignment operator.

  /// Assigns \a rhs to \a *this.  If dimensions
  /// agree, we reuse the data vector, otherwise we construct a new
  /// matrix object and swap ourselves with it.
  /// \todo Efficiency improvements for special cases.

  template< typename tag >
  matrix< T , dense_column > const& operator=( matrix< T , tag > const& rhs ) ;

private:

  /// Data vector.

  std::vector< T > data ;

} ;


/// Specialization for dense matrices.

template< typename T >
struct matrix< T , dense > {

  typedef T      &       reference ;
  typedef T const& const_reference ;

  matrix( long const rows , long const columns = 1 )
  : rows_   ( ( assert( rows    >= 1 ) , rows    ) ) ,
    columns_( ( assert( columns >= 1 ) , columns ) ) ,
    data( rows * columns )
  {}

  /// Default constructor makes an empty matrix (zero rows and columns).

  matrix()
  : rows_( 0 ) , columns_( 0 ) , data( 0 ) {}

  /// Templated copy constructor delegating to assignment operator.
  /// \todo This is somewhat inefficient.  The members are initialized and
  /// immediately overwritten.

  template< typename tag >
  matrix( matrix< T , tag > const& A )
  : rows_( 0 ) , columns_( 0 ) , data( 0 )
  { *this = A ; }

  /// Swap with other matrix

  void swap( matrix< T , dense >& ) ;

  /// Set all elements to specified value.

  void fill( T const& ) ;

  long rows   () const { return rows_    ; }
  long columns() const { return columns_ ; }

  reference operator()( long i , long j = 1 )
  { check_indices( *this , i , j ) ;
    return data[ row_major_index( i , j , columns() ) ] ; }

  const_reference operator()( long i , long j = 1 ) const
  { check_indices( *this , i , j ) ;
    return data[ row_major_index( i , j , columns() ) ] ; }


  /// Templated assignment operator.

  /// Assigns \a rhs to \a *this.  If dimensions
  /// agree, we reuse the data vector, otherwise we construct a new
  /// matrix object and swap ourselves with it.
  /// \todo Efficiency improvements for special cases.

  template< typename tag >
  matrix< T , dense > const& operator=( matrix< T , tag > const& rhs ) ;

private:

  /// Construct a column vector from \a x.

  matrix( std::vector< T > const& x )
  : rows_( x.size() ) , columns_( 1 ) , data( x ) {}

  long rows_    ;
  long columns_ ;

  /// Data vector.

  std::vector< T > data ;

} ;



/// Specialization for fixed matrices.  Dimensions are fixed at compile
/// time.
///
/// Note that this type does NOT support efficient swapping!

template< typename T , long rows_ , long columns_ >
struct matrix< T , fixed< rows_ , columns_ > > {

  typedef T      &       reference ;
  typedef T const& const_reference ;

  /// Default constructor sets all elements to zero.

  matrix()
  { assert( rows_ >= 1 ) ; assert( columns_ >= 1 ) ; fill( T() ) ; }

  void fill( T const& t ) 
  { std::fill( data , data + rows_ * columns_ , t ) ; }

  /// Construct from an initializer object.

  template< typename R >
  matrix( initializer< R > const& init )
  { assert( init.size == rows_ * columns_ ) ;
    std::copy( init.data , init.data + init.size , data ) ; }

  /// Templated copy constructor delegating to assignment operator.

  template< typename tag >
  matrix( matrix< T , tag > const& A )
  { *this = A ; }

  long rows   () const { return rows_    ; }
  long columns() const { return columns_ ; }

  reference operator()( long i , long j = 1 )
  { check_indices( *this , i , j ) ;
    return data[ row_major_index( i , j , columns() ) ] ; }

  const_reference operator()( long i , long j = 1 ) const
  { check_indices( *this , i , j ) ;
    return data[ row_major_index( i , j , columns() ) ] ; }


  /// Templated assignment operator.

  /// Assigns \a rhs to \a *this.  Dimensions must agree.

  template< typename tag >
  matrix< T , fixed< rows_ , columns_ > > const&
  operator=( matrix< T , tag > const& rhs )
  { return assign( *this , rhs ) ; }

private:

  /// Data array.

  T data[ rows_ * columns_ ] ;

} ;


/// Specialization for transposed matrix.

template< typename T , typename tag >
struct matrix< T , transposed< tag > > {

  /// Type of contained matrix.
  typedef matrix< T , tag > contained_type ;

  /// Constant element reference.
  typedef T const& const_reference ;

  /// Constructor: initialize contained matrix.
  matrix( contained_type const& A ) : A( A ) { }

  long rows   () const { return A.columns() ; }
  long columns() const { return A.rows   () ; }

  const_reference operator()( long i , long j = 1 ) const
  { return A( j , i ) ; }

  contained_type const& contained() const { return A ; }

private:

  /// Reference to contained matrix.
  contained_type const& A ;

} ;


/// View of a std::vector< T > as a matrix.

template< typename T >
struct matrix< T , vector > {

  matrix( long const rows , long const columns , std::vector< T > const& v )
  : rows_( rows ) , columns_( columns ) , v( v )
  { assert( static_cast< unsigned long >( rows * columns ) == v.size() ) ; }

  long rows   () const { return rows_    ; }
  long columns() const { return columns_ ; }

  T const& operator()( long i , long j = 1 ) const
  { check_indices( *this , i , j ) ;
    return v[ row_major_index( i , j , columns() ) ] ; }

private:

  long const rows_    ;
  long const columns_ ;

  std::vector< T > const& v ;

} ;


/// View of a std::vector< std::vector< T > > as a matrix.

template< typename T >
struct matrix< T , vector_vector > {

  /// Initialize and assert that all rows have the same size.
  matrix( std::vector< std::vector< T > > const& A ) ;

  long rows() const { return rows_ ; }
  long columns() const { return columns_ ; }

  T const& operator()( long i , long j = 1 ) const
  { check_indices( *this , i , j ) ; return A[ i - 1 ][ j - 1 ] ; }

private:

  /// Reference to contained matrix.
  std::vector< std::vector< T > > const& A ;

  long rows_ ;  ///< Number of rows.
  long columns_ ; ///< Number of columns.

} ;



/// Specialization representing a matrix with a unary operator
/// applied to each element.

/// The bandwidths are the same as for the contained matrix.  Note that
/// this requires that f(x) = 0 whenever x=0 and f is defined.  In
/// particular, this holds if f(x)=ax for some a or f(x)=1/x.

template< typename T , typename tag , typename op >
struct matrix< T , un_op< tag , op > > {

  /// Type of contained matrix.
  typedef matrix< T , tag > contained_type ;

  /// Constructor initializing contained matrix and operator.
  matrix( contained_type const& A , op const& f ) : A( A ) , f( f ) {}

  /// \return Number of rows.
  long rows() const { return A.rows() ; }
  /// \return Number of columns.
  long columns() const { return A.columns() ; }

  const T operator()( long i , long j = 1 ) const
  { return f( A( i , j ) ) ; }
  
  contained_type const& operand () const { return A ; }
  op             const& function() const { return f ; }

private:

  /// The contained matrix.
  contained_type const& A ;
  /// The operator, contained by value.
  const op f ;

} ;



///
/// A matrix defined by a function of its indices.  
/// 

template< typename T , typename op , long rows_ , long columns_ >
struct matrix< T , index_fun< op , rows_ , columns_ > > {

  // f must define lower_bandwidth() and upper_bandwidth() functions!
  // NOTE f is stored by value.

  matrix( op const& f ) : f( f ) {}

  long rows()    const { return rows_    ; }
  long columns() const { return columns_ ; }

  T const operator()( long i , long j = 1 ) const
  { check_indices( *this , i , j ) ; return f( i , j ) ; }

  op const& function() const { return f ; }

private:

  op const f ;

} ;


/// Specialization representing an immutable submatrix 
/// defined by a range of rows.

/// This expression template provides a const operator() for accessing
/// the contained matrix.

template< typename T , typename tag >
struct matrix< T , const_row_range< tag > > {

  /// Type of the contained matrix.

  typedef matrix< T , tag > contained_type ;


  /// Row submatrix of A from row index first to last (incl.)

  matrix( contained_type const& A , long first , long last )
  : A( A ) ,
    shift( first - 1 ) ,
    rows_( last - shift )
  { assert( 1 <= first && first <= last && last <= A.rows() ) ; }

  /// \return Number of rows.
  long rows() const { return rows_ ; }
  /// \return Number of columns.
  long columns() const { return A.columns() ; }

  /// \return The row index in the contained matrix, given the row index
  /// in the submatrix.
  long original_row( long i ) const
  { assert( 1 <= i && i <= rows() ) ; return i + shift ; }

  T const& operator()( long i , long j = 1 ) const
  { return A( original_row( i ) , j ) ; }

private:

  /// The contained matrix.
  contained_type const& A ;

  /// Row shift value.

  /// If i is a row index in the contained matrix, and ii the
  /// corresponding row index in the submatrix, we have ii = i + shift.
  /// \sa original_row.

  const long shift ;

  /// Number of rows.

  const long rows_ ;

} ;


/// Specialization representing a mutable submatrix defined by a range of rows.

/// This expression template provides both const and non-const operator()
/// for accessing the contained matrix.  It also provides an assignment
/// operator for assigning to the row range.

template< typename T , typename tag >
struct matrix< T , row_range< tag > > {

  /// Type of the contained matrix.

  typedef matrix< T , tag > contained_type ;


  /// Row submatrix of A from row index first to last (incl.)

  matrix( contained_type& A , long first , long last )
  : A( A ) ,
    shift( first - 1 ) ,
    rows_( last - shift )
  { assert( 1 <= first && first <= last && last <= A.rows() ) ; }

  /// \return Number of rows.

  long rows() const { return rows_ ; }

  /// \return Number of columns.

  long columns() const { return A.columns() ; }

  /// \return The row index in the contained matrix, given the row index
  /// in the submatrix.

  long original_row( long i ) const
  { assert( 1 <= i && i <= rows() ) ; return i + shift ; }

  T& operator()( long i , long j = 1 )
  { return A( original_row( i ) , j ) ; }

  T const& operator()( long i , long j = 1 ) const
  { return A( original_row( i ) , j ) ; }

  /// Assign right hand side to *this.

  template< typename tag2 >
  matrix< T , row_range< tag > > const&
  operator=( matrix< T , tag2 > const& rhs ) ;

  /// Set all elements to specified value.

  void fill( T const& ) ;

private:

  /// The contained matrix.

  contained_type& A ;

  /// Row shift value.

  /// If i is a row index in the contained matrix, and ii the
  /// corresponding row index in the submatrix, we have ii = i + shift.
  /// \sa original_row.

  const long shift ;

  /// Number of rows.

  const long rows_ ;

} ;


/// Specialization for a row scaling.

/// This object represents the matrix diag( S ) * A.  S must be a column
/// vector (i.e., columns() = 1).

template< typename T , typename tag1 , typename tag2 >
struct matrix< T , row_scaled< tag1 , tag2 > > {

  /// Type of the scaling vector S.
  typedef matrix< T , tag1 > contained_type_1 ;
  /// Type of the scaled matrix A.
  typedef matrix< T , tag2 > contained_type_2 ;

  /// Initialize const references.
  matrix( contained_type_1 const& S , contained_type_2 const& A )
  : S( S ) , A( A )
  { assert( S.rows() == A.rows() && S.columns() == 1 ) ; }

  /// \return Number of rows.
  long rows() const { return A.rows() ; }
  /// \return Number of columns.
  long columns() const { return A.columns() ; }

  T operator()( long i , long j = 1 ) const
  { return S( i , 1 ) * A( i , j ) ; }

private:

  /// The scaling vector.
  contained_type_1 const& S ;
  /// The scaled matrix.
  contained_type_2 const& A ;

} ;



/// Specialization representing result of elementwise binary operation.

/// Typically, f will be plus or minus.  We assert that f(0,0)=0 so that
/// the bandwidths of the result are bounded by the bandwidth of the
/// operands.

template< typename T , typename tag1 , typename tag2 , typename op >
struct matrix< T , bin_op< tag1 , tag2 , op > > {

  /// Type of the first operand.
  typedef matrix< T , tag1 > contained_type_1 ;
  /// Type of the second operand.
  typedef matrix< T , tag2 > contained_type_2 ;

  /// Initialize \a A, \a B, and \a f and assert that \a f(0,0)=0 and dimensions agree.
  matrix( contained_type_1 const& A , contained_type_2 const& B , op const& f )
  : A( A ) , B( B ) , f( f )
  { assert( A.columns() == B.columns() && A.rows() == B.rows() ) ;
    assert( f( 0 , 0 ) == 0 ) ; }

  /// \return Number of rows.
  long rows() const { return A.rows() ; }
  /// \return Number of columns.
  long columns() const { return A.columns() ; }

  T operator()( long i , long j = 1 ) const
  { return f( A( i , j ) , B( i , j ) ) ; }

  contained_type_1 const& left_operand () const { return A ; }
  contained_type_2 const& right_operand() const { return B ; }

private:

  contained_type_1 const& A ;
  contained_type_2 const& B ;

  /// Binary operation, NOTE: stored by value.
  op f ;

} ;


/// Specialization representing the product of two matrices.
template< typename T , typename tag1 , typename tag2 >
struct matrix< T , product< tag1 , tag2 > > {

  typedef matrix< T , tag1 > contained_type_1 ;
  typedef matrix< T , tag2 > contained_type_2 ;

  /// Initialize with operands A and B.
  matrix( contained_type_1 const& A , contained_type_2 const& B )
  : A( A ) , B( B )
  { assert( A.columns() == B.rows() ) ; }

  long rows()    const { return A.rows   () ; }
  long columns() const { return B.columns() ; }

  T operator()( long i , long j = 1 ) const ;

  contained_type_1 const& left_operand () const { return A ; }
  contained_type_2 const& right_operand() const { return B ; }

private:

  contained_type_1 const& A ;
  contained_type_2 const& B ;

} ;

//@}


/// \defgroup range_funcs Functions defining nonzero row and column ranges.
///
//@{

/// First nonzero column for given row.

template< typename T , typename tag >
inline long first_column( matrix< T , tag > const& A , long i ) {

  cpl::math::check_row_index( A , i ) ;
  return std::max( i - lower_bandwidth( A ) , 1L ) ;

}

/// Last nonzero column for given row.

template< typename T , typename tag >
inline long last_column( matrix< T , tag > const& A , long i ) {

  cpl::math::check_row_index( A , i ) ;
  return std::min( i + upper_bandwidth( A ) , A.columns() ) ;

}


/// First nonzero row for given column.

template< typename T , typename tag >
inline long first_row( matrix< T , tag > const& A , long j ) {

  cpl::math::check_column_index( A , j ) ;
  return std::max( j - upper_bandwidth( A ) , 1L ) ;

}

/// Last nonzero row for given column.

template< typename T , typename tag >
inline long last_row( matrix< T , tag > const& A , long j ) {

  cpl::math::check_column_index( A , j ) ;
  return std::min( j + lower_bandwidth( A ) , A.rows() ) ;

}

//@}


/// \defgroup proxy_factories Functions returning matrices.
///
//@{

/// \return The column vector constructed from x.

template< typename T >
inline matrix< T , dense > column( std::vector< T > const& x )
{ return matrix< T , dense >( x ) ; }


/// \return The transposed of A.

template< typename T , typename tag >
inline const matrix< T , transposed< tag > >
transpose( matrix< T , tag > const& A )
{ return matrix< T , transposed< tag > >( A ) ; }


/// \return An assignable row submatrix.

/// \param A The original matrix.
/// \param first The first row index in A.
/// \param last The last row index in A.

template< typename T , typename tag >
inline matrix< T , row_range< tag > >
rows( matrix< T , tag >& A , long first , long last )
{ return matrix< T , row_range< tag > >( A , first , last ) ; }

/// \return A constant row submatrix.

/// \param A The original matrix.
/// \param first The first row index in A.
/// \param last The last row index in A.

template< typename T , typename tag >
inline const matrix< T , const_row_range< tag > >
rows( matrix< T , tag > const& A , long first , long last )
{ return matrix< T , const_row_range< tag > >( A , first , last ) ; }


/// \return The row scaling diag(S)*A.

template< typename T , typename tag1 , typename tag2 >
inline const matrix< T , row_scaled< tag1 , tag2 > >
row_scaling( matrix< T , tag1 > const& S , matrix< T , tag2 > const& A )
{ return matrix< T , row_scaled< tag1 , tag2 > >( S , A ) ; }

/// \return The matrix ( f( A( i , j ) ) ).

template< typename T , typename tag , typename op >
inline const matrix< T , un_op< tag , op > >
apply( matrix< T , tag > const& A , op const& f )
{ return matrix< T , un_op< tag , op > >( A , f ) ; }

/// \return -A (unary minus).

template< typename T , typename tag >
inline const matrix< T , un_op< tag , std::negate< T > > >
operator-( matrix< T , tag > const& A )
{ return apply( A , std::negate< T >() ) ; }

/// \return \a aA (scalar multiplication).

template< typename T , typename tag >
inline const
matrix< T , un_op< tag , std::binder1st< std::multiplies< T > > > >
operator*( T const& a , matrix< T , tag > const& A )
{ return apply( A , std::bind1st( std::multiplies< T >() , a ) ) ; }


/// \return \a AB (matrix multiplication).

template< typename T , typename tag1 , typename tag2 >
inline const matrix< T , product< tag1 , tag2 > >
operator*( matrix< T , tag1 > const& A , matrix< T , tag2 > const& B )
{ return matrix< T , product< tag1 , tag2 > >( A , B ) ; }


/// Helper for identity matrix
struct kronecker_delta {

  int operator()( long const i , long const j = 1 ) const
  { return i == j ; }

  long upper_bandwidth() const { return 0 ; }
  long lower_bandwidth() const { return 0 ; }

} ;

/// \return \a rows * \a columns identity matrix.

template< typename T , long rows , long columns >
inline const matrix< T , index_fun< kronecker_delta , rows , columns > >
identity()
{ return matrix< T , index_fun< kronecker_delta , rows , columns > >(
    kronecker_delta() ) ; }


// Element-wise compare of two matrices.
template< typename T , typename tag1 , typename tag2 >
bool const
operator==( matrix< T , tag1 > const& A , matrix< T , tag2 > const& B ) { 
  
  if( A.rows   () != B.rows   () ) { return false ; } 
  if( A.columns() != B.columns() ) { return false ; }

  for( long i = 1 ; i <= A.rows   () ; ++i )
  for( long j = 1 ; j <= A.columns() ; ++j ) 
  { if( A( i , j ) != B( i , j ) ) { return false ; } }

  return true ;

}


template< typename T , typename tag1 , typename tag2 >
bool const
operator!=( matrix< T , tag1 > const& A , matrix< T , tag2 > const& B ) 
{ return !( A == B ) ; }

  

/// \return A+B (matrix addition).

template< typename T , typename tag1 , typename tag2 >
inline const matrix< T , bin_op< tag1 , tag2 , std::plus< T > > >
operator+( matrix< T , tag1 > const& A , matrix< T , tag2 > const& B )
{ return
  matrix< T , bin_op< tag1 , tag2 , std::plus< T > > >
  ( A , B , std::plus< T >() ) ; }


/// \return ( A( i , j ) * B( i , j ) ) (element-wise multiplication)

template< typename T , typename tag1 , typename tag2 >
inline const matrix< T , bin_op< tag1 , tag2 , std::multiplies< T > > >
elementwise_multiply(
  matrix< T , tag1 > const& A ,
  matrix< T , tag2 > const& B
)
{ return
  matrix< T , bin_op< tag1 , tag2 , std::multiplies< T > > >
  ( A , B , std::multiplies< T >() ) ; }


/// \return A - B (matrix subtraction).

template< typename T , typename tag1 , typename tag2 >
inline const matrix< T , bin_op< tag1 , tag2 , std::minus< T > > >
operator-( matrix< T , tag1 > const& A , matrix< T , tag2 > const& B )
{ return
  matrix< T , bin_op< tag1 , tag2 , std::minus< T > > >
  ( A , B , std::minus< T >() ) ; }


template< typename T , typename tag1 , typename tag2 >
inline const matrix< T , tag1 >&
operator+=( matrix< T , tag1 >& A , matrix< T , tag2 > const& B ) 
{ return A = A + B ; }

template< typename T , typename tag >
inline matrix< T , tag > const& 
operator*=( matrix< T , tag >& A , T const& s ) 
{ return A = s * A ; }

template< typename T , typename tag1 , typename tag2 >
inline const matrix< T , tag1 >&
operator-=( matrix< T , tag1 >& A , matrix< T , tag2 > const& B ) 
{ return A = A - B ; }


//@}

/// \defgroup scalar_return Functions taking matrix arguments and returning scalars.
///
//@{

/// \return max | A( i , j ) |

template< typename T , typename tag >
T max_abs_norm( matrix< T , tag > const& A ) ;

/// \return sum | A( i , j ) |

template< typename T , typename tag >
T sum_abs_norm( matrix< T , tag > const& A ) ;

/// \return sum A( i , j )^2

template< typename T , typename tag >
T square_sum( matrix< T , tag > const& A ) ;

/// \return sqrt( sum A( i , j )^2 )

template< typename T , typename tag >
T norm_2    ( matrix< T , tag > const& A )
{ return std::sqrt( square_sum( A ) ) ; }

/// \return Minimum element of A

template< typename T , typename tag >
T min( matrix< T , tag > const& A ) ;

/// \return Maximum element of A

template< typename T , typename tag >
T max( matrix< T , tag > const& A ) ;


/// \return Inner product of x and y.

/// \param x A column vector.
/// \param y A column vector of the same dimension.

template< typename T , typename tag1 , typename tag2 >
T const 
inner_product( matrix< T , tag1 > const& x , matrix< T , tag2 > const& y ) ;

//
// An alias for the inner product.  Always use with parentheses: ( x | y )
//

template< typename T , typename tag1 , typename tag2 >
inline T const
operator|( matrix< T , tag1 > const& x , matrix< T , tag2 > const& y ) 
{ return inner_product( x , y ) ; }


// 
// The cross-product of two three-dimensional column vectors.
//

template< typename T , typename tag1 , typename tag2 >
matrix< T , fixed< 3 , 1 > > const
cross_product( 
  matrix< T , tag1 > const& x ,
  matrix< T , tag2 > const& y
) {

  assert( x.rows() == 3 ) ;
  assert( y.rows() == 3 ) ;
  
  assert( x.columns() == 1 ) ;
  assert( y.columns() == 1 ) ;

  matrix< T , fixed< 3 , 1 > > z ;

  z( 1 , 1 ) = x( 2 , 1 ) * y( 3 , 1 ) - x( 3 , 1 ) * y( 2 , 1 ) ;
  z( 2 , 1 ) = x( 3 , 1 ) * y( 1 , 1 ) - x( 1 , 1 ) * y( 3 , 1 ) ;
  z( 3 , 1 ) = x( 1 , 1 ) * y( 2 , 1 ) - x( 2 , 1 ) * y( 1 , 1 ) ;

  return z ;

}


//
// Determinant of a 3x3 matrix.
//
// This is generated code from a little Mathematica script.
//

template< typename T , typename tag >
T const det_3( matrix< T , tag > const& A ) {
  
  assert( 3 == A.rows   () ) ;
  assert( 3 == A.columns() ) ;

  return
    -A(1,3)*A(2,2)*A(3,1) + A(1,2)*A(2,3)*A(3,1) + A(1,3)*A(2,1)*A(3,2) -
     A(1,1)*A(2,3)*A(3,2) - A(1,2)*A(2,1)*A(3,3) + A(1,1)*A(2,2)*A(3,3) ;

}


//
// Inverse of a 3x3 matrix.
//
// This is generated code from a little Mathematica script.
//

template< typename T , typename tag >
matrix< T , tag > const
inverse_3( matrix< T , tag > const& A ) {

  assert( A.rows   () == 3 ) ;
  assert( A.columns() == 3 ) ;

  T const D = det_3( A ) ;

  assert( D != 0 ) ;

  matrix< T , fixed< 3 , 3 > > B ;

  B( 1 , 1 ) = 1 / D * ( -A(2,3)*A(3,2) + A(2,2)*A(3,3) ) ;
  B( 1 , 2 ) = 1 / D * (  A(1,3)*A(3,2) - A(1,2)*A(3,3) ) ;
  B( 1 , 3 ) = 1 / D * ( -A(1,3)*A(2,2) + A(1,2)*A(2,3) ) ;
  B( 2 , 1 ) = 1 / D * (  A(2,3)*A(3,1) - A(2,1)*A(3,3) ) ;
  B( 2 , 2 ) = 1 / D * ( -A(1,3)*A(3,1) + A(1,1)*A(3,3) ) ;
  B( 2 , 3 ) = 1 / D * (  A(1,3)*A(2,1) - A(1,1)*A(2,3) ) ;
  B( 3 , 1 ) = 1 / D * ( -A(2,2)*A(3,1) + A(2,1)*A(3,2) ) ;
  B( 3 , 2 ) = 1 / D * (  A(1,2)*A(3,1) - A(1,1)*A(3,2) ) ;
  B( 3 , 3 ) = 1 / D * ( -A(1,2)*A(2,1) + A(1,1)*A(2,2) ) ;

  return B ;

}


//@}



/// \defgroup misc Utility functions.
///
//@{

/// Swap A and B by calling swap member function.

/// Note that this only works for the basic matrix types.

template< typename T , typename tag >
void swap( matrix< T , tag >& A , matrix< T , tag >& B ) ;

/// Assign delegate.  Set \a LHS to \a RHS.

/// No redimensioning of the LHS is done, we assert() that the
/// dimensions agree.  To be called by the assignment operators.
/// We iterate across the nonzero band of \a LHS.
/// \return \a LHS.

template< typename T , typename tag1 , typename tag2 >
matrix< T , tag1 > const&
assign( matrix< T , tag1 >& LHS , matrix< T , tag2 > const& RHS ) ;



//
// Convenience conversion functions.
//

// Return an initializer containing i elements.
// Precondition:  a is a T[] of i elements.

template< typename T >
inline initializer< T >
initialize( T const& a )
{ return initializer< T >( a , sizeof( a ) / sizeof( a[ 0 ] ) ) ; }


/// Column vector initializer.
///
/// NOTE.  T is derived from argument types, you need to use e.g.
/// doubles in order to initialize a 3-vector of doubles!

template< typename T >
inline matrix< T , fixed< 3 , 1 > > const
column_vector( T const& v1 , T const& v2 , T const& v3 ) {

  matrix< T , fixed< 3 , 1 > > ret ;
  ret( 1 , 1 ) = v1 ;
  ret( 2 , 1 ) = v2 ;
  ret( 3 , 1 ) = v3 ;
  
  return ret ;

}


/// Return a matrix proxy using row major indexing into a vector.

template< typename T >
inline matrix< T , vector > as_matrix
( long const rows , long const columns , std::vector< T > const& v )
{ return matrix< T , vector >( rows , columns , v ) ; }


/// Return a matrix proxy from a vector of vectors of T

template< typename T >
inline matrix< T , vector_vector > as_matrix
( std::vector< std::vector < T > > const& A ) 
{ return matrix< T , vector_vector >( A ) ; }


/// Apply f to each element in the nonzero band of A.

template< typename T , typename tag , typename op >
void apply_nonzero( matrix< T , tag > const& A , op& f ) ;


/// Convert a column vector to a std::vector.

template< typename T , typename tag >
std::vector< T > to_vector( matrix< T , tag > const& c ) ;


} // end namespace math

} // end namespace cpl


//
// Template definitions.
//

inline unsigned long cpl::math::row_major_index
( long const i , long const j , long const columns ) {

  assert( i >= 1       ) ;
  assert( j >= 1       ) ;
  assert( j <= columns ) ;

  return static_cast< unsigned long >( ( i - 1 ) * columns + j - 1 ) ;

}


template< typename T , typename tag , typename op >
void cpl::math::apply_nonzero
( cpl::math::matrix< T , tag > const& A , op& f ) {

  assert( A.rows   () >= 1 ) ;
  assert( A.columns() >= 1 ) ;

  for( long i = 1 ; i <= A.rows() ; ++i ) {

    for( long j  = cpl::math::first_column( A , i ) ;
              j <= cpl::math::last_column ( A , i ) ;
            ++j )
    { f( A( i , j ) ) ; }

  }

}



namespace {


template< typename T >
struct sum_abs_accumulator {

  sum_abs_accumulator() : res( T() ) {}

  void operator()( T const& t ) { res += std::fabs( t ) ; }

  T result() const { return res ; }

private:

  T res ;

} ;

template< typename T >
struct max_abs_accumulator {

  max_abs_accumulator() : res( T() ) {}

  void operator()( T const& t )
  { if( std::fabs( t ) > res ) { res = std::fabs( t ) ; } }

  T result() const { return res ; }

private:

  T res ;

} ;


template< typename T >
struct square_sum_accumulator {

  square_sum_accumulator() : res( T() ) {}

  void operator()( T const& t )
  { res += t * t ; }

  T result() const { return res ; }

private:

  T res ;

} ;


template< typename T >
struct max_accumulator {

  max_accumulator() : first( true ) {}

  void operator()( T const& t )
  { if( first ) { res = t ; first = false ; }
    else        { if( t > res ) res = t ; } }

  T result() const { assert( !first ) ; return res ; }

private:

  T res ;
  bool first ;

} ;


template< typename T >
struct min_accumulator {

  min_accumulator() : first( true ) {}

  void operator()( T const& t )
  { if( first ) { res = t ; first = false ; }
    else        { if( t < res ) res = t ; } }

  T result() const { assert( !first ) ; return res ; }

private:

  T res ;
  bool first ;

} ;


} // anonymous namespace


template< typename T , typename tag >
T cpl::math::sum_abs_norm( cpl::math::matrix< T , tag > const& A ) {

  sum_abs_accumulator< T > saa ;
  cpl::math::apply_nonzero( A , saa ) ;
  return saa.result() ;

}

template< typename T , typename tag >
T cpl::math::square_sum( cpl::math::matrix< T , tag > const& A ) {

  square_sum_accumulator< T > a ;
  cpl::math::apply_nonzero( A , a ) ;
  return a.result() ;

}

template< typename T , typename tag >
T cpl::math::max_abs_norm( cpl::math::matrix< T , tag > const& A ) {

  max_abs_accumulator< T > maa ;
  cpl::math::apply_nonzero( A , maa ) ;
  return maa.result() ;

}

template< typename T , typename tag >
T cpl::math::max( cpl::math::matrix< T , tag > const& A ) {

  max_accumulator< T > ma ;
  cpl::math::apply_nonzero( A , ma ) ;
  if( cpl::math::has_out_of_band_entries( A ) ) { ma( T() ) ; }
  return ma.result() ;

}

template< typename T , typename tag >
T cpl::math::min( cpl::math::matrix< T , tag > const& A ) {

  min_accumulator< T > ma ;
  cpl::math::apply_nonzero( A , ma ) ;
  if( has_out_of_band_entries( A ) ) { ma( T() ) ; }
  return ma.result() ;

}


template< typename T >
void cpl::math::matrix< T , cpl::math::dense >::swap
( cpl::math::matrix< T , dense >& other ) { 
  std::swap( rows_ , other.rows_ ) ;
  std::swap( columns_ , other.columns_ ) ;
  std::swap( data , other.data ) ;

}


template< typename T , typename tag >
inline void cpl::math::swap( 
  cpl::math::matrix< T , tag >& A , 
  cpl::math::matrix< T , tag >& B 
) 
{ A.swap( B ) ; }


template< typename T >
void cpl::math::matrix< T , cpl::math::dense >::fill( T const& x ) {

  assert( rows() >= 1 ) ;
  assert( columns() >= 1 ) ;

  std::fill( data.begin() , data.end() , x ) ;

}


template< typename T >
cpl::math::matrix< T , cpl::math::vector_vector >
::matrix( std::vector< std::vector< T > > const& A )
: A( A ) {

  rows_ = A.size() ;

  if( rows_ == 0 ) { columns_ = 0 ; return ; }

  assert( rows_ > 0 ) ;

  columns_ = A[ 0 ].size() ;

  assert( columns_ > 0 ) ;

  for( long i = 2 ; i <= rows_ ; ++i )
  { assert( static_cast< long >( A[ i - 1 ].size() ) == columns_ ) ; }

}


template< typename T , typename tag1 , typename tag2 >
T const cpl::math::inner_product(
  cpl::math::matrix< T , tag1 > const& x ,
  cpl::math::matrix< T , tag2 > const& y
) {

  long n = x.rows() ;
  assert( n == y.rows() ) ;
  assert( x.columns() == 1 ) ;
  assert( y.columns() == 1 ) ;

  T ret = T() ;

  for( long i = 1 ; i <= n ; ++i ) {

    ret += x( i , 1 ) * y( i , 1 ) ;

  }

  return ret ;

}


template< typename T , typename tag1 , typename tag2 >
T cpl::math::matrix< T , cpl::math::product< tag1 , tag2 > >
::operator()( long i , long j ) const {

  using namespace cpl::math ;

  check_indices( *this , i , j ) ;

  T ret = T() ;

  for(
    long k =  std::max( first_column( A , i ) , first_row( B , j ) ) ;
         k <= std::min( last_column( A , i ) , last_row( B , j ) ) ;
       ++k )

    ret += A( i , k ) * B( k , j ) ;

  return ret ;

}


template< typename T , typename tag1 , typename tag2 >
const cpl::math::matrix< T , tag1 >&
cpl::math::assign(
  cpl::math::matrix< T , tag1 >      & A ,
  cpl::math::matrix< T , tag2 > const& B
) {

  assert( A.rows() == B.rows() ) ;
  assert( A.columns() == B.columns() ) ;

  for( long i = 1 ; i <= A.rows() ; ++i ) {

    for( long j =  cpl::math::first_column( A , i ) ;
              j <= cpl::math::last_column ( A , i ) ;
            ++j ) {

      A( i , j ) = B( i , j ) ;

    }

  }

  return A ;

}


template< typename T >
template< typename tag >
cpl::math::matrix< T , cpl::math::dense > const&
cpl::math::matrix< T , cpl::math::dense >::operator=
( cpl::math::matrix< T , tag > const& rhs ) {

  if(    rhs.rows()    != rows()
      || rhs.columns() != columns() ) {

    cpl::math::matrix< T , cpl::math::dense >
      ret( rhs.rows() , rhs.columns() ) ;

    cpl::math::assign( ret , rhs ) ;
    this->swap( ret ) ;

  } else { cpl::math::assign( *this , rhs ) ; }

  return *this ;

}


template< typename T >
template< typename tag >
inline
cpl::math::matrix< T , cpl::math::dense_column > const&
cpl::math::matrix< T , cpl::math::dense_column >::operator=
( cpl::math::matrix< T , tag > const& rhs ) {

  assert( 1 == rhs.columns() ) ;
  this->data.resize( rhs.rows() ) ;
  cpl::math::assign( *this , rhs ) ;

}


template< typename T , typename tag >
template< typename tag2 >
cpl::math::matrix< T , cpl::math::row_range< tag > > const&
cpl::math::matrix< T , cpl::math::row_range< tag > >::operator=
( cpl::math::matrix< T , tag2 > const& rhs ) {

  assert( rows() == rhs.rows() ) ;
  assert( columns() == rhs.columns() ) ;

  assign( *this , rhs ) ;

  return *this ;

}


template< typename T , typename tag >
void cpl::math::matrix< T , cpl::math::row_range< tag > >::fill
( T const& x ) {

  assert( rows() >= 1 ) ;
  assert( columns() >= 1 ) ;

  for( long i = 1 ; i <= rows() ; ++i ) {

    for( long j = 1 ; j <= columns() ; ++j ) {

      ( *this )( i , j ) = x ;

    }

  }

}


template< typename T , typename tag >
std::vector< T > cpl::math::to_vector
( cpl::math::matrix< T , tag > const& c ) {

  assert( c.columns() == 1 ) ;

  std::vector< T > ret( c.rows() ) ;
  for( unsigned long i = 0 ; i < ret.size() ; ++i )
  { ret[ i ] = c( i + 1 , 1 ) ; }

  return ret ;

}


#endif // CPP_LIB_MATRIX_H
